/*
  File:    LeastSquares.cpp
  Purpose: Utilities to find the least squares best fit polynomial for a set 
           of data points. 

  Contact: Paul Macklin
           pmacklin@math.uci.edu
		  http://math.uci.edu/~pmacklin
*/

#include <cstdio>
#include <cstdlib>
#include <iostream>

using namespace std;

#include "../ImprovedMath.h"
#include "../Matrix.h"
#include "Polynomial.h"

Polynomial LeastSquaresPolynomialFit( Matrix& Xpoints, Matrix& Ypoints, int Degree )
{
 Polynomial Output( Degree );
 int m = Xpoints.TellCols();
 if( Ypoints.TellCols() < m || m < Degree+1 )
 {
  cout << "Least squares error: Not enough data points given!" << endl; 
  return Output;
 }
 
 Matrix Coefficients(Degree+1,Degree+1);
 Matrix RHS(1,Degree+1);
 
 int i,j,k;
 for( j=0 ; j <= Degree ; j++ )
 {
  for( i=0 ; i < m ; i++ )
  { *RHS(j) += ( *Ypoints(i) )*IntPow( *Xpoints(i) , j); }
   
  for( k=0 ; k <= Degree ; k++ )
  {
   for( i=0 ; i < m ; i++ )
   { *Coefficients(j,k) += IntPow(*Xpoints(i),j+k); }
  }
 }
 
 Matrix temp = GaussianEliminationPartialPivoting( Coefficients, RHS );

 for( i=0; i <= Degree ; i++ )
 { Output.SetCoefficient( i , *temp(i) ); }
 return Output;
}
